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ABSTRACT 
In 1951, G.W. Brown proposed an iterative algorithm called ‘fictitious play’ 


for solving two-person zero-sum games. Although it is an effective method, the 


fictitious play algorithm converges slowly to the value of the game. Recently, Gass, . 


Zafra, and Qiu proposed a modification that applies to symmetric games, i.e., games 
with skew-symmetric payoff matrices. To solve non-symmetric games via their 
modification, the games must be made symmetric via a transformation. 

Gass, Zafra, and Qiu reported that their modified algorithm converges faster 
than the original fictitious play on a collection of randomly generated games. 
However, their results on non-symmetric games only apply to games whose values 
are near zero. When game values are far away from zero, this thesis empirically 
shows that the original fictitious play algorithm can outperform the modified one. 

Gass, Zafra, and Qiu’s method is static, in that the symmetric transformation 
is done once prior to the start of their modified algorithm. However, they suggested 
the exploration of dynamic methods where the transformation is periodically 
revised. This thesis proposes and investigates the convergence behavior of one 


dynamic transformation technique for solving general two-person zero-sum games. 
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I. INTRODUCTION | 

When it was proposed in 1951, Brown’s fictitious play algorithm was not known 
to converge, but it was applied to a few two-person zero-sum (TPZS) games with 
success. Robinson [Ref. 1] later proved its convergence. Despite this result, the 
algorithm is not always the method of choice for solving TPZS games. In their book, 
Szép and Forg6 [Ref. 2] state that ‘Computational experience available up to now 
indicates that, for the solution of general matrix games, linear programming is the most 
efficient method.’ However, the fictitious play algorithm is not impractical. For some 
large games or games where it is impossible or impractical to enumerate all possible 
strategies a priori (See, e.g., Ref. [3]), the fictitious play algorithm may be the only 
effective choice. 

Recently, Gass, Zafra, and Qiu (GZQ) [Ref. 4] proposed a modification to the 
fictitious play algorithm. The modification assumes that the game is symmetric. Since 
non-symmetric games can be transformed into symmetric ones (see, e.g., Ref. [5]), their 
modification also applies to non-symmetric games. In their paper, GZQ also report that 
their modification converges faster than the original fictitious play algorithm. Their 
results are based on random games with 100x100 payoff matrices whose elements are 
Uniform random numbers between —100 and 100. 


Table 1.1: Values of Games with Uniform[-100,100] Payoffs 


Density Game values among a sample of 50 games 
|_minimum | average _|_ maximum _ 











Table 1.1 summarizes the values of some games with Uniform [-100,100] 
payoffs. In the table, the size of the payoff matrices is 100x100, i.e., each player has 100 
strategies available, and the density is varied from 20% to 100%. For each combination 
of size and density, 50 random games are generated and solved by the simplex method 
(see, e.g., Ref. [6]) to find the game values. The minimum, average, and maximum game 
values for each sample of 50 games are reported in Table 1.1. They indicate that the 
values of the games with Uniform [-100,100] payoffs are near zero. Thus, GZQ tests are 
mainly on games with values that are near zero. 

One objective of this thesis is to consider a larger class of random games and 
numerically compare the original against the modified fictitious play algorithm as 
proposed by GZQ. In addition, GZQ’s modification is static, in that the symmetric 
transformation is performed once prior to start of their algorithm. This thesis also 
proposes an alternate modification in which the transformation is periodically updated 
with a different scaling parameter. The numerical results reported here indicate that this 
new modification converges more rapidly than both the original and modified fictitious 
play procedures proposed by GZQ. 

To make the thesis self-contained, Chapter II describes symmetric and non- 
symmetric TPZS games. Chapter III discusses the existing techniques for solving these 
games and proposes an alternate modification for the fictitious play algorithm. In 
Chapter IV, variations of the fictitious play algorithm are compared against the original. 


Finally, Chapter V summarizes the results. 











II. TWO-PERSON ZERO-SUM GAMES 

A. DEFINITION 

A two-person zero-sum (TPZS) game is a situation where there are two decision- 
makers (players) having directly opposite interests. Ina TPZS game, one player, P1, has 
m Strategies available and the other, P2, has n strategies. The outcome of the game 
depends on the strategy used by each player. If P1 and P2 choose the in and j" strategies, 
respectively, then the outcome of the game, denoted as a,;, represents the amount P2 has 
to pay P1. In other words, the payoffs to P1 and P2 are aj and -aj, respectively. Note that 
the sum of the two payoffs is zero, which explains the name of the game. 

A TPZS game is completely defined when the payoff for each pair of P1 and P2 
strategies is determined. These payoff can be summarized as an mxn matrix, generally 


referred to as a ‘payoff matrix’, i.e., 


Qa, 42 Qin 
A a1 O22 a2, 
ani an? — Qin 


In playing the game, both players are assumed to choose the best strategy to achieve the 
most favorable outcome. This means that P1, the player receiving the payoff a, would 
choose the strategy that maximizes the outcome of the game. On the other hand, P2, 
having to pay aj, would choose the strategy that minimizes the outcome instead. 

When a TPZS game has an equilibrium point, each player can guarantee an — 
outcome of the game by always choosing one particular strategy. When equilibrium can 


not be achieved, both players are assumed to randomize their strategies. In other words, 





P1 and P2 choose strategies i and j, respectively, independently with probability x; and y,. 
The randomized strategies, x = (x, ..., Xm) and y = (yj, ..., Yn)’ must satisfy the 
following conditions: 

1) O<x, <1 for alli=1,....m 


ii) O<y <1 for allj=1,...,n 
jit) Six] 
iv) y,=1 
To choose the best randomized strategy, P1 must solve the following optimization 
max { min). ¥,4,]; yx =1, x,20,i=1,....m} (2.1) 
j=l i=! 
Similarly, P2 must solve the following problem to obtain the best randomized strategy 


min { max [ )/a,y,]; y=, y,20,j21,....7) (2.2) 
j=l j 


Let x* and y* denote the optimal strategies for P1 and P2. Then, v” =(x")’ Ay’ is the 


value of the game. In addition, it is well-known (see, e.g., Ref. [6]) that 


v= max { min [ )/ x2, ] Dae: =1, x,20,i=1,...,m}, and 


i=] i=] 
v’ = min {max[)/a,yj}) Doyj=1, yj 20,j=1,-.0} 
Jel j 


B. SYMMETRIC GAMES 
A TPZS game is said to be symmetric if its payoff matrix is skew-symmetric, 1.e., 


A=-A' or aj = -a; for all i andj. The value of a symmetric game is always zero (see 








[Ref.7]). In addition, the optimal randomized strategies for P1 and P2 are the same, i.e., 
x =y. 
C. NON-SYMMETRIC GAMES 

When the payoff matrix is not skew-symmetric, the game is non-symmetric and 
the value of the game may be nonzero. However, it is possible to transform a non- 
symmetric game into a symmetric one. In the literature, there are two transformation 
techniques, one proposed by von Neuman and the other by Gale, Kuhn and Tucker 
(GKT) [Ref. 5]. Among these two techniques, the latter is more attractive 
computationally, for the resulting symmetric game has a much smaller payoff matrix. 


Given a payoff matrix A, the GKT technique defines the following payoff matrix 


0 A -¢, 
S=|-A’™ 0  ¢, 
it ann 


Where c; and c2 are vectors in R” and R", respectively. Furthermore, the components of 
both c; and c2 are the same and they all are equal to 6 , a positive constant. The 0’s in $ 
represent zero matrices of compatible dimension. 

When the original payoff matrix A is of size mxn, the resulting oayoff matrix S is 
of size (m+n+1)x (m+n+1). Clearly, S is skew-symmetric. Let (p*,g",/°) be an optimal 


randomized strategy for the symmetric game with payoff matrix S. If the original game 


has a positive value, it is possible to show (see [Ref. 5]) that )) p; = >, q; and the 


i=] j=l 


optimal strategies for P1 and P2 in the original game are x" =~ p” and y’ =+-q° where 


= wre . In addition, the value of the original game is 6/°/ nu. 


j=l 








For games with zero or negative values, it is possible to ‘scale’ or add a constant 
w to the payoff matrix A and obtain A(w) = [a; + w]. When w is sufficiently large, the 


resulting game has a positive value. 








I. SOLVING TWO-PERSON ZERO-SUM GAMES 


A. LINEAR PROGRAMMING 


The problem for determining an optimal randomized strategy for P1 in equation 


2.1 is equivalent to the following linear program (see, e.g., Ref. [6]): 


OPT-1: maximize vy 


subject to i a,x, 2v forjel,....n 
i 


x, 20 fori=1,....m 
Similarly, the problem for determining an optimal randomized strategy for P2 in equation 
2.2 1s equivalent to the following linear program: 
OPT-2: igi w 
subject to >a); <w fori=1,....m 
j 
j 
yy 20 forj=1,...,n 
It is easy to show that problems OPT-1 and OPT-2 are dual of each other. Moreover, if 
(v*, x*) and (w*, y*) are optimal to problems OPT-1 and OPT-2, respectively, then v* = 
w*, 
B. REGULAR FICTITIOUS PLAY 
To adopt the convention used in GZQ, the fictitious play algorithm as proposed 
by G. W. Brown in Ref. [8] is henceforth referred to as the regular fictitious play 
algorithm or RFP. To describe RFP, let A; and A’ denote the i” row and the j" column of 


the payoff matrix A. As before, x = (x1. ..., Xj, ..., Xm) and y = (yj, ..., Yj, --+5 Yay 





represent the randomized strategies for Pl and P2, respectively. After k fictitious plays, 


define 


1) Uk) =([U\(k), ...,Ui(k),..., U,Ak)\ as the cumulative payoff vector for P1, 


li) aj(k) as the number of times P1 uses strategy i in k plays, 


ili) V(k) = [V(A), ...,-Vj(K), «5 VARY" as the cumulative payoff vector for P2, 


iv) bj(k) as the number of times P2 uses strategy j in k plays, 


v) m(k) as the lower bound for the value of the game, and 


vi) M(k) as the upper bound for the value of the game. 


Below, RFP is stated with a stopping tolerance of ¢ > 0. 


Step 0: 


Step 1: 


Step 2: 


Step 3: 


Regular Fictitious Play Algorithm (RFP) 
Set U(0) = 0, V(0) = 0, a(0) = 0 for all i=1,...,m, b(O) =0 for all j = 1, ...,, 
m(Q) =~ 0, M(O) = 0, and k= 1. Go to Step 1. 
Let r = arg max; {U,(k-1),..., Ui(kK-1),..., Un(k-1)} (break ties arbitrarily). Set 
Vk)’ = V(k-1)' + A, and ak) = afk — 1) + 1. 1fk > 1, compute M(K) = min{M(k 
— 1), g@nU,(k —-1) } and go to Step 2. 
Let s = arg min;{V;(k),..., Vj(k), .... VaCk)} (break ties arbitrarily). Set U(k) = 
U(k-1) + A’ and b,(k) = b,(k — 1) + 1. Compute m(k) = max{m(k - 1), £V,(k) } 
and go to Step 3. 
If [M(k) — m(k)] < &, stop. The randomized strategy pair, x = 
1(a,(k),..., a, (k))” and y = 1(b,(k),..., b, (k))" , 1S approximately optimal. 


Otherwise, set k =k + 1 and go to Step 1. 


In the first occurence Step 1, the r strategy (i.e, r" row of the payoff matrix A) is 


arbitrarily assigned to P1 since U(Q) is a zero vector. In Step 2, P2 minimizes the 











cumulative payoff to P1 by choosing the s* strategy. If the gap, i.e., the difference 
between the upper and lower bounds, in Step 3 is sufficiently small, the algorithm 
computes the resulting randomized strategies for both players and terminates. 


To illustrate, consider the payoff matrix 


4 6 8 
A= , 
ree 


Table 3.1 summarizes the first ten iterations of RFP. In Step 1, P1 arbitrarily chooses 
Strategy 2 (or row 2) since V() = 0 and V(1) is set to (7, 5, 1)’. In Step 2, P2 then 
chooses strategy 3 (or column 3) since, against V(1), P2 only has to pay 1 unit to P1 using 
this strategy. This sets U(1) = (8, 1)’. The sequence of play just described is 
summarized in the first row (k = 1) of Table 3.1. The remaining rows are similarly 
obtained. In the fourth row (k = 4), m(4) is the same as m(3) since it is larger that V,(4)/4 
= 4.75. Similarly, in the ninth row (k = 9), M(9) = M(8) since M(8) is less than U>(k- 
1)/(k-1) = 44/8 = 5.5 


Table 3.1: Ten Iterations of the Regular Fictitious Play Algorithm 


SE ee m{(k) 
0 0 7 wae 






qOomMmAINKNRWHN 
fmd eed pre here ped pet tet SC) 1) 
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Cc. MODIFIED FICTITIOUS PLAY 
There are three variations of RFP discussed here. The first variation is called the 


modified fictitious play or MFP algorithm. This algorithm only applies to symmetric 








games. The other two are intended for non-symmetric games. Both use the GKT 
technique discussed in Chapter II to transform a non-symmetric game into a symmetric 
one. The variation described in GZQ, referred to here as ‘Static MFP’ or SMFP, 
performs the symmetric transforms once prior to the start of the algorithm. The 
remaining variation, ‘Dynamic MFP’ or DMEFP, is new, Although, the basic idea is 
alluded to in GZQ. In this new variation, the transformation is periodically updated with 
new lower bounds on the value of the original game. As demonstrated in the next 
chapter, this variation converges faster on a collection of randomly generated games. 
1. Modified Fictitious Play for Symmetric Games 
The MFP algorithm for symmetric games is essentially the same as RFP. There is 
an additional step in MFP to take advantage of the fact that the value of a symmetric 
game is always zero and the optimal randomized strategies for P1 and P2 are the same. 
For an nxn payoff matrix, the MFP algorithm requires a switching interval K in addition 
to the stopping tolerance ¢ and it can be stated as follows. 
Modified Fictitious Play Algorithm (MFP) 
Step 0: Set U(0) = 0, V(O) = 0, a,(0) = 0 for all 7= 1, ...,n, B(O) =0 for all 7 = 1, 
....2, m(O) =— 0, M(O) =o, andk= 1. Goto Step 1. 
Step 1: Letr=arg max; {U;(k-1),..., U(k-1),..., U,(k—-1)} (break ties arbitrarily). 
Set V(k)’ = V(k-1)’ + A, and a,{k) =a,(k—1)+1. Ifk> 1, then set M(k) = 


min{M(k — 1),—+-U_(k —1) }. Go to Step 2. 


> (k-1) 
Step 2: Lets =arg minj{V,(A), ..., V(x), .... Vi(k)} (break ties arbitrarily). Set 
U(k) = U(k-1) + A’, bs(k) = b(k — 1) + 1, and m(k) = max{m(k — 1), 


+V,(k) }. Go to Step 3. 


10 





Step 3: = If [M(k) — m(k)] < €, stop and either the randomized strategy x 
=1(a,(k),...,a,(k))7 ory = 4(b,(k),....b,(k))* 18 approximately optimal. 
Otherwise, go to Step 4. 

Step 4: If mod(k, K)>0, setk=k-+1 and goto Step 1. Otherwise, set B= 


min{Imin{LV,(k)}1, Imax{4U,(k)}1}. If 8 =|min{LV,(k)}|, then set 


U(k) =- va) and b(k) = a(k). Otherwise, set V(k) = — U(k) and a(k) = 
b(k). Setk=k+1 and go to Step 1. 
Steps 0 to 3 are essentially the same as those in RFP. The notation used in MFP reflects 
the fact that the payoff matrix, A, is square. In Step 4, MFP computes at every K 


iterations two game value estimates, Imin{;V, (k)}| and Imax{4U, (k)}|. Depending on 
J i 


which estimate is closer to zero, the cumulative payoff vectors are ‘switched’ and both 
a(k) and b(k) are made equal to reflect the better of the two game value estimates. As 
empirically demonstrated in GZQ, the best value for K is one. 

Zi Modified Fictitious Play with Static Scaling 

For a non-symmetric game, Static MFP or SMFP first transforms the mxn payoff 


matrix into the following skew-symmetric (m + n + 1) x (m +n + 1) payoff matrix: 


0 A(w) -c¢, 
Sw)=|A(w)’ 0 c, 


ci -c, (0 
where A(w) = [a,;; + w] for some scaling factor w > 0 and the components of c, € R” and 


and c2 € R’ are all equal to 5 > 0. For the above transformation to be valid, the scaling 


factor w must be large enough to ensure that the payoff matrix, A(w), yields a game with 





a positive value. In GZQ, the scaling parameter w is set to Imina;/+1. Then, SMFP is 


essentially MFP applied to S(w). 


As before, let S,(w) and S°(w) denote the r™ row and the s“ column of the payoff 


matrix S(w). Then, the Static MFP algorithm can be stated as follows. 


Step 0: 


Step 1: 


Step 2: 


Step 3: 


Modified Fictitious Play Algorithm with Static Scaling 
(SMFP) 


Set U(0) = 0, V(0) = 0, a(0) = 0 for all i= 1, ..., (mtn+1), b(O) =0 for all j = 1, 
., (m+n+1), m(O) = — 0, M(0) = 0, and k = 1. Go to Step 1. 

Let r = arg max; {U,(kK-1),..., Ui(kK-1),..., U¢mansiy(K-1)} (break ties arbitrarily). 

Set V(k)’ = V(k-1)" + S{w) and a,(k) = a(k — 1) + 1. Go to Step 2. 

Let s = arg min,{ V(A), ..., Vj(k), .--, Vinen+1)(k)} (break ties arbitrarily). Set U(k) 

= U(k-1) + S°(w) and b,(k) = b(k — 1) + 1. Go to Step 3. 


Set 


x, = a fori=1,..., m, 
75a 


y; — for j = 1,...,n, 
aes (k) 


for i= 1,..., m, and 


for j= 1,...,7. 


Then, compute 


m(k) = max{m(k —1), min > A, (w)x;,min D A,(w) p,} and 


j=l 4 


M (k) = min{M (k -1), max ) A; (w)y;,max D7, A,(w)q ;} 





If [M(k) — m(&)] < €, stop and either the randomized strategy pair (x, y) or (p, q) 
is approximately optimal to the game with payoff matrix A. Otherwise, go to 
Step 4. 

Step 4: If mod(k, K) > 0, setk=k+1 and goto Step 1. Otherwise, set f= 


min{Imin{;V,(k)}1, lmax{;U,(k)}1}. If # =!min{;V,(k)}I, then set U(k) =— 


V(k) and b(k) = a(k). Otherwise, set V(k) = — U(k) and a(k) = b(k). Setk=k+1 
and go to Step 1. 

The main difference between MFP and SMFP is in Step 3 where the upper and 
lower bounds for the value of the game are calculated. The bounds in Step 3 of SMFP 
are for the original game, not the one with S(w). 

Step 3 may involve dividing by zero, but MATLAB, the numeric computational 
software in which the algorithm is implemented, handles this eventuality correctly (see 
appendix E). 

3. Modified Fictitious Play with Dyaanle Scaling 

The results in Table 3.2 show that setting w to lmina,|l + 1 makes SMFP converge 
slowly. In fact, a better value for w is —v, where v is the value of the game. Table 3.2 
displays the results on eight random games with 25%, 50%, 75%, and 100% density. 
Four games have 10x 10 payoff matrices and the other have 20x20. Elements of all 
payoff matrices are uniformly distributed between —100 and 100. The stopping tolerance, 
€, is 0.1 and the switching interval, K,1s 1. The value of the game is obtained by solving 


the corresponding linear program discussed earlier. 
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Table 3.2: Numerical Results for MFP with two scaling parameters 


Matrix Matrix | Game w = Imin(q;,)l + 1 
Size Density | value CPU | Iteration CPU Iteration 
(sec) (sec) 


| 10x10 | 100% _| -14.73 | 12.36] 3496 | 6.81 [2092 | 














In Table 3.2, the CPU times for MFP with w = — v are between 10% to 80% of 
those for MFP with w = Imina,l+1. Over all eight games, there is a 44% reduction in 
CPU times when w =- v. A similar conclusion also holds for the number of iterations. 

The above results motivate the idea of adjusting the scaling parameter 
periodically. To be valid, the resulting payoff matrix A(w) must yield a positive game 
value. One method is to simply set w to the negative of the best lower bound, i.e., 

w =—m(k). 

Let L be the rescaling interval. Then, the Dynamic MFP algorithm or DMFP can 

be stated as follows: 


Modified Fictitious Play Algorithm with Dynamic Scaling 
(DMFP) 


Steps 0 - 3: Same as SMFP 


Step 4: If mod(k, K) > 0, go to Step 5. Otherwise, set B= min{Imin{;V,(k)}l, 
J 
Imax{1U,(k)}l}. If B = Imin{;-V,(k)}I, then set U(k) = — V(k) and b(k) = 
i J 


a(k). Otherwise, set V(k) =— U(k) and a(k) = b(k) and go to Step 5. 
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Step5: Setk=k+1. Ifmod(k, L) >0, goto Step 1. Otherwise, set w = — m(k), 


recompute S(w), and go to Step 1. 





IV. NUMERICAL RESULTS 

A. DATA GENERATION 

To compare the algorithms discussed in Chapter III, random games with 100x 100 
payoff matrices are generated. For symmetric games, the payoffs, aj, are uniformly 
distributed between —100 and 100 (see appendix A). For non-symmetric games, three 
groups of payoff matrices are considered. In Group 1, the nonzero payoffs are uniformly 
distributed between —100 and 100. For Group 2, they are between —200 and 0. Finally, 
aj Ss for Group 3 are between 0 and 200. 

To generate payoff matrices with density 0, where 0 < @ < 1, a Uniform random 
number, j1;;, between O and 1 is generated for each pair of strategies (7, j). If py < @, then 


a Uniform random number is generated for aj. Otherwise, a;;= 0. The games in Groups 


1, 2, and 3 differ by a constant if and only if 0 = 1. 


B. PARAMETER VALUES 

All algorithms are terminated when the gap is less than or equal to 0.1, 1-e., the 
stopping tolerance, ¢, equals 0.1. As suggested by GZQ, the switching interval K is 1 for 
MFP, SMEFP, and DMFP. 

For SMFP, 6 is set to 0.75(max(a;)-min(a;)). Recall that the parameter 6 is the 
constant used in defining c; and c2 for the symmetric transformation. In our preliminary 
study, other values for 6, e.g., 0.25(max(a;)-min(a,)), 0.50(max(a;;)-min(a,j)), (max(ay)- 
min(ajj)), ai 1, did not show significant improvement in CPU time. 

For games in Groups 1 and 2, the scaling parameter w is Imin(a;)| + 1. Since 


games in Group 3 already have positive values, w is simply set to zero. 
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The parameter settings for DMFP are the same as those for SMFP. The additional 
parameter for the dynamic version is the rescaling interval, L, which is set at 5000. 
C. IMPLEMENTATION 

All four algorithms were implemented as MATLAB (Version 5.0) functions. The 
programs for all four functions and game generation are listed in the appendices. The 
reported CPU times in the following sections were generated using a 300 MHz Pentium 
II personal computer with 64 MEG of RAM. 
D. SYMMETRIC GAMES 

The results in Table 4.1 show that MFP outperforms RFP on four symmetric 


games from Group 1 with various densities. 


Table 4.1: Computational Results for Symmetric Games in Group 1 


Method Gap Lower Upper CPU Iterations 
bound bound (sec) (xX 10°) 
1 
] 


0.1000 | -0.0502 | 0.0498 40.37 | 0.0080 
0.1000 | -0.0500 | 0.0499 63.55 | 0.0102 
0.0998 | -0.0505 | 0.0493 83.37 | 0.0179 
0.1000 | -0.0532 | 0.0468 
0.0994 | -0.0497 | 0.0497 79.92 | 0.0161 


The CPU times for MFP in Table 4.1 are between 7.5% and 10.0% of those for RFP. 





Similarly, MFP also uses fewer iterations; they are between 4% and 5% of those for RFP. 
Figure 4.1 displays a typical convergence behavior of RFP and MFP on 
symmetric games. MFP seems to converge directly to the solution. On the other hand, 
RFP has a long convergence tail. The success of MFP can be attributed to the switching 
of cumulative payoffs and strategies in Step 4, which takes full advantage of the special 


solution structure of symmetric games. 
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Figure 4.1: Convergence Behavior of RFP and MFP for a Symmetric Game 
in Group 1 with 100% Density 


E. NON-SYMMETRIC GAMES 

Two sets of experiments were performed, one set to compare RFP against SMFP 
and the other to compare RFP against DMFP. In these experiments, all three algorithms 
are terminated as soon as a 0.1 gap is achieved. As demonstrated below, DMFP is 
superior to RFP and RFP is superior to SMFP. 

1, RFP and SMFP 

Tables 4.2 to 4.4 summarize results for games in Groups 1, 2, and 3, respectively. 
Each table reports results on four random games, each with different densities. With the 
exception of game 3 in Table 4.4, RFP takes between 20% to 75% less time than SMFP 
to achieve a gap of 0.1 or better. For game 3 in Table 4.4, the CPU time of RFP is 


approximately 3% more than that of SMFP. In general, SMFP requires fewer iterations. 
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However, one iteration of SMFP is much more computationally intensive than one 
iteration of RFP. 


Table 4.2: Computational Results for Games in Group 1 
pT ee Ss ee 
bound bound (sec) (x 10°) 
SMFP 0.0997 -0.7103 | -0.6106 1466.79 | 0.2402 
SMFP 0.0999 -0.2627 | -0.1628 1617.66 | 0.2788 
SMFP 0.0999 0.5136 0.6135 2247.34 | 0.2898 
SMFP 0.0998 0.1217 0.2216 1795.85 | 0.2486 

Table 4.3: Computational Results for Games in Group 2 
ape eer Tes Tee ee 
bound bound (sec.) (x 10°) 
[=| ge | ee Sees ae | 
SMFP 0.0997 | -24.8027 | -24.7029 | 2849.15 | 0.5316 
SMFP 0.0999 | -49.6847 | -49.5847 | 2183.56 | 0.4056 
SMFP 0.0999 | -73.8156 | -73.7157 | 1578.72 | 0.3191 
0.1000 -97.6296 | 798.40| 0.5140 
SMFP 0.0998 | -97.7225 | -97.6226 | 1398.89 | 0.2880 

Table 4.4: Computational Results for Games in Group 3 
SS 
bound bound (sec) (x 10°) 
aed 
SMFP 0.0999 25.8965 | 25.9964 2850.68 | 0.4885 
=e 
SMFP 0.0995 48.6532 | 48.7528 1269.44 | 0.2354 
: 
SMFP 0.0998 73.4502 | 73.5500 1133.66 | 0.2399 


100% RFP 0.0999 | 99.0131 | 99.1130 825.47 | 0.5287 
SMFP 0.0998 | 99.0190 | 99.1188 | 1568.51 | 0.2869 























It is also interesting to note in Tables 4.3 and 4.4 that SMFP takes the most CPU 
time in reaching a 0.1 gap for the game with 25% dense payoff matrices. 
Figures 4.2 to 4.4 graphically display a typical convergence behavior for RFP and 


SMEFP for games in Groups 1, 2 and 3 when ¢ (the stopping tolerance) = 0.1. These 
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figures show that RFP converges to the game value faster than SMFP in terms of CPU 


time for all three groups of games. 
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Figure 4.2: Convergence Behavior of RFP and SMFP for a Game in Group 1 with 50% Density 
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Figure 4.3: Convergence Behavior of RFP and SMFP for a Game in Group 2 with 50% Density 
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Figure 4.4: Convergence Behavior of RFP and SMFP for a Game in Group 3 with 50% Density 


2. RFP and DMFP 

As in the above subsection, Tables 4.5 to 4.7 summarize the computational results 
for games in Groups 1, 2, and 3, respectively. Each table reports results on four random 
games, each with different densities. In all 12 games, DMFP takes between 36% to 67% 
less time than RFP to achieve a gap of 0.1 or better. As with the static version, DMFP 


requires fewer iterations than RFP. However, unlike its counterpart, small numbers of 


DMF? iterations do not translate into more CPU seconds. 








Table 4.5: Computational Results for Games in Group 1 


i a a= 
bound bound (sec) (x 10°) 

alg) ee 
(DMP | 0.0999] -0-7177 | -0.6178 | 159.45] 0.0228 

50% |__RFP | 0.1000 | -0.8660 | -0.7660 | 461.92] 0.2959 
75% |__RFP | 0.1000 | 03411 | 02411 | 564.47] 0.3576 
joo% |__RFP__| 0.1000 | 0.1049 | 0.2048 | 617.20] 03892 











Table 4.6: Computational Results for Games in Group 2 


ap = [ee ee os Pee 
bound bound (sec) (x 10°) 

25% | ___RFP | 0.10 | -23.89 | -23.79 | 906.16] 05449 _ 
50% |___RFP__| 010 | -#9.42 | 49.32 | 1029.64) 0.6765_ 
75% |__RFP___|_0.10_| -7159 | -71.49 | 103688] 05775 
100% | __RFP__| 0.10 | -100.02 | 99.92 | 657.13] 0.4192_ 
[_DMrP [0.10 | -100.02 | -99.93 [378.05] 0.0652_ 















Table 4.7: Computational Results for Games in Group 3 


Tee ee ee 
bound bound (sec) (x 10°) 

5% |__RFP | 010 | 23.05 | 23.15 | 80048 05308 
[_DMFP__[0.10__[23.05_|-23.15_| 263.59] 0.0401 

ld gg ee ee 
Ge gd ee ea 
















0.1 

0.1 

DMFP 

[00% | __RFP___| 010 _| 9780 | 9790 | 111186 
pMFP__| 0.10 | 97.80 | 9790 | 642.08 


Figures 4.5 to 4.7 graphically illustrate typical RFP and DMFP convergence 
behavior for games in Groups 1, 2 and 3, respectively. As before, ¢, the stopping 
tolerance, is set at 0.1. As the numerical results suggest, DMFP outperforms RFP which, 


in turn, outperforms SMFP. 
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Figure 4.5: Convergence Behavior of RFP and DMFP 
for a Game in Group 1 with 50% Density 
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Figure 4.6: Convergence Behavior of RFP and DMFP 
for a Game in Group 2 with 50% Density 
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Figure 4.7: Convergence Behavior of RFP and DMFP 
for a Game in Group 3 with 50% Density 
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Vv. CONCLUSIONS AND SUGGESTION FOR FURTHER WORK 

This thesis proposes an alternate modification to the fictitious play algorithm 
which can be considered as a dynamic variation of the one proposed by GZQ. The 
modification proposed here is dynamic, in that the symmetric transformation is updated 
periodically with a different scaling parameter. The results in Chapter IV indicate that 
this dynamic variation is computationally advantageous when compared to the original 
fictitious play or the static variation proposed by GZQ. 

The results in Chapter IV also confirm that GZQ’s modified fictitious play 
algorithm outperforms the original when applied to symmetric games. The results are 
reversed for non-symmetric games. The original fictitious play algorithm is rather robust 
and outperforms GZQ’s static modification on random non-symmetric games. 

Finally, this thesis also identifies several topics for further investigation. One 
topic is to investigate the choice of 6 in defining c, and cz in the sated matrix S(w). 
The other is to investigate other choices isi w, the scaling parameter, perhaps in 


conjunction with the choice of 6. 
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APPENDIX A. MATLAB CODES FOR GENERATING SKEW-SYMMETRIC 
MATRICES FOR SYMMETRIC GAMES 


function A = skewSym(m,d) 
% This function is used for generating an (m Xm) skew-symmetric matrix 
% with density d percent. — 
A = zeros(m,m); 
for i= 1:m 
for j = 1:m 
if i<j 
if d > round(rand(1)*100) % Generate elements of A ~ U(-100,100) 
if 0.5 > rand(1) 
A(i,j) = ceil(rand(1)*100); 
else 
A(i,j) = floor(-100*rand(1)); 
end 
end 
end 
ifi>j 
A(i,j) sa -A(j,1); 
end 
end 
end 
return 
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_ APPENDIX B. MATLAB CODES FOR GENERATING MATRICES FOR NON- 
SYMMETRIC GAMES 


function A = density(m,n,d) 
% This function generates an (m Xn) general matrix with d % density 
A = zeros(m,n); % A is a payoff matrix 
for i= 1:m 
for } = I:n 
if d > round(rand(1)*100) % Generate elements of A ~ U(-100,100) 
if 0.5 > rand(1) 
A(i,j) = ceil(rand(1)*100); 
else 
A(i,j) = floor(-100*rand(1)); 
end 
end 
end 
end 
return 
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APPENDIX C. MATLAB CODES FOR RFP ALGORITHM 


function [Upb,Lwb,gp,gap,k,compTime,reqflops,P1,P2] = RFP(gapneed,A) 
% This function performs the RFP and gives the result whenever it achieves the 
% needed gap 


[m,n] = size(A); 

V = zeros(1,n); % V is a P2's (choices)cumulative payoff vector 
U = zeros(m,1); % U is a P1's (choices)cumulative payoff vector 
Xx = zeros(m,1); % x is a P1's strategy vector 

y = zeros(1,n); % y is a P2's strategy vector 

r=0; 

mx = 10000; 


t = cputime; 
f = flops; 


% The first iteration (k = 1) 

first = ceil(rand(1)*m); % first is the i th row that P1 arbitrarily chooses 
V=V + A(first,:); 

x(first) = x(first) + 1; 

lower(1) = min(V); 


second = find(V==min(V)); % second is the j th column that P2 chooses to respond 

if length(second) > 1 
nn = ceil(rand(1)*(length(second)); 
second = second(nn); 

end 

U =U +A(:,second); 

y(second) = y(second) + 1; 

upper(1) = max(U); 

gap(1) = upper(1) - lower(1); 

k=1; 

ko = 1; 


while gap(ko) > gapneed 
k=k+1; 
ko = ko + 1; 
st = find(U==max(U)); % st is the i th row giving max payoff that P1 expect 
if length(st) > 1 
nn = ceil(rand(1)*length(st)); 
st = st(nn); 
end 
V=V+ ACst,:); 
x(st) = x(st) + 1; 
second = find(V==min(V)); % second is the j th column that P2 chooses to respond 
if length(second) > 1 
nn = ceil(rand(1)*length(second)); 
second = second(nn); 
end 
U=U + A(,second); 
y(second) = y(second) + 1; 
upper(ko) = max(U)/k; 
lower(ko) = min(V)/k; 
if upper(ko) > upper(ko-1) 
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upper(ko) = upper(ko-1); 

end 

if lower(ko) < lower(ko-1) 
lower(ko) = lower(ko-1); 

end 

gap(ko) = upper(ko) - lower(ko); 


% In order to save the memory space, 
% the following loop controls the result vectors (lower, upper and gap) 
% to have at most mx = 10,000 elements inside. 
if mod(k-1,mx) == 0 
r=r+1; 
disp({r, gap(ko-1)}) %display the gap at every 10,000 iterations 
lower(1) = lower(ko); 
lower = lower(1:ko-1); 
upper(1) = upper(ko); 
upper = upper(1:ko-1); 
gap(1) = gap(ko); 
gap = gap(1:ko-1); 
ko = 1; 
end 
en 
k; 
if k > mx 
if ko < 5000 
lower = [lower((mx-5000+ko+1):mx),lower(1:ko)]; 
upper = [upper((mx-5000+ko+1):mx),upper(1:ko)]; 
gap = [gap((mx-5000+ko+1):mx),gap(1:ko)]; 
ko = 5000; 
else 
lower = lower(1:ko); 
upper = upper(1:ko); 
gap = gap(1:ko); 
end 
else 
lower = lower(1:ko); 
upper = upper(1:ko); 
gap = gap(1:ko); 
end 
Pl = x./k; 
P2 = y./k; 
Upb = upper(ko); Lwb = lower(ko); 
gp = gap(ko); 
reqflops = flops - f; 
compTime = cputime - t; 


return 
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APPENDIX D. MATLAB CODES FOR MFP ALGORITHM FOR SYMMETRIC 
GAMES 


function [k,Upb,Lwb,gp,gap,compTime,reqflops] = MFPSym(gapneed,B) 

% This function performs the MFP with an attempt to duplicate MFP algorithm in Gass, 
% Zafra's paper. 

% This function gives the result whenever it achieves the needed gap 


K=1; 

r= 0; mx = 10000; 

A =B; % input matrix game 

[m,n] = size(A); 

V = zeros(1,n); % V is a P2's (choices)cumulative payoff vector 
U = zeros(m,1); % U 1s a PI's (choices)cumulative payoff vector 
x = zeros(m,1); % x 1s a PI's strategy vector 

y = zeros(1,n); % y is a P2's strategy vector 

t = cputime; 

f = flops; 


% The first iteration (k = 1) 

first = round(1+rand(1)*(m-1)); % first is the i th row that P1 arbitrarily chooses 
V=V+A(fiirst,:); 

x(first) = x(first) + 1; 

lower(1) = min(V); 


second = find(V==min(V)); % second is the j th column that P2 chooses to respond 
if length(second) > 1 
nn = round(1+rand(1)*(length(second)-1)); 
second = second(nn); 
end 
U =U + AC,,second); 
y(second) = y(second) + 1; 
upper(1) = max(U); 
gap(1) = upper(1) - lower(1); 


k= I: 
ko = 1; 
if K == 
if abs(max(U)) < abs(min(V)) 
V=-(U); 
else 
U=-(V); 
end 
end 
while gap(ko) > gapneed 
k=k+1; 
ko = ko +1; 


st = find(U==max(U)); % st is the i th row giving max payoff that P1 expect 
if length(st) > 1 
nn = round(1+rand(1)*(length(st)-1)); 
st = st(nn); 
end 
V=V+ A(st,:); 
x(st) = x(st) + 1; 
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second = find(V==min(V)); % second is the j th column that P2 chooses to respond 
if length(second) > 1 
nn = round(1+rand(1)*(length(second)-1)); 
second = second(nn); 
end 
U =U + A(:,second); 
y(second) = y(second) + 1; 
upper(ko) = max(U)/k; 
Jower(ko) = min(V)/k; 
if upper(ko) > upper(ko-1) 
upper(ko) = upper(ko-1); 
end 
if lower(ko) < lower(ko-1) 
lower(ko) = lower(ko-1); 
end 
gap(ko) = upper(ko) - lower(ko); 


if mod(k,K) == 0 
if abs(max(U)/k) < abs(min(V)/k) 


V=-(U); 
else 
U=-(V); 
end 
end 


% In order to save the memory space, 
% the following loop controls the result vectors (lower, upper and gap) 
% to have at most mx = 10,000 elements inside. 
if mod(k-1,mx) == 
r=r+1; 
disp([r, gap(ko-1)]}) %display the gap at every 10,000 iterations 
lower(1) = lower(ko); 
lower = lower(1:mx); 
upper(1) = upper(ko); 
upper = upper(1:mx); 
gap(1) = gap(ko); 
gap = gap(1:mx); 
ko= 1; 
end 
end 


reqflops = flops - f; 
compTime = cputime - t; 


if k > mx 
if ko < 5000 %number of iterations 1s in [r*10000 , r*10000+5000) 
lower = [lower((mx-5000+ko+1):mx),lower(1:ko)]; 
upper = [upper((mx-5000+ko+1):mx),upper(1:ko)]J; 
gap = [gap((mx-5000+ko+1):mx),gap(1:ko)]; 
ko = 5000; 
else %number of iterations is in [r*10000+5000 , (r+1)*1 0000) 
lower = lower(1:ko); 
upper = upper(1:ko); 
gap = gap(1:ko); 
end 
else %number of iterations < 10000 
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lower = lower(1:ko); 
upper = upper(1:ko); 
gap = gap(1:ko); 

end 7 

Upb = upper(end); Lwb = lower(end); gp = gap(end); 

return 
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APPENDIX E. MATLAB CODES FOR SMFP 


function [count, Upb,Lwb,gp,gap,compTime,reqflops,P1,P2] = SMFP(gapneed,A,w) 

% This function performs the MFP with an attempt to duplicate MFP algorithm in Gass, 
% Zafra's paper. 

% This function gives the result whenever it achieves the needed gap 


r= 0; 

mx = 10000; 

k= 1: 

[m,n] = size(A); 

alpha = 0.75; 

c = alpha*(max(max(A)) - min(min(A))); 
An=A+w; 


S = [zeros(m) An -c*ones(m,1); 
-An’ zeros(n) c*ones(n,1); . 
c*ones(1,m) -c*ones(1,n) 0]; 
[sm,sn] = size(S); 


U = zeros(sm, 1); 
V = zeros(1,sn); 
X = zeros(sm, 1); 
y = zeros(1,sn); 


t = cputime; 
f = flops; 


%start the first iteration 

st = ceil(rand(1)*sm); % Row player randomly choose his row strategy 
x(st) = x(st) + 1; | 

V=V+S(st:); 


nd = find(V == min(V)); 

if (length(nd) > 1) % Randomly choose to break the tie 
i = ceil(rand(1)*length(nd)); 
nd = nd(i); 

-end 

y(nd) = y(nd) + 1; 

U=U + S(,nd); 


Jo To To To To To To Fo Yo To To No Yo To To To To To Yo Yo 


xbarl = x(1:m)./sum(x(1:m)); % This part will cause zero division errors 
ybarl = x(m+1:m+n)./sum(x(m+1l:m+n)); % for some initial iterations, however MATLAB 
xbar2 = y(1:m)./sum(y(1:m)); % can continue to the end of mission 


ybar2 = y(m+1 :m+n)./sum(y(m+1:m+n)); 


lower(1) = max(min(xbarl'*A), min(xbar2*A)); 
upper(1) = min(max(A*ybarl), max(A*ybar2’)); 
gap(1) = upper(1) - lower(1); 


ik==1 
if (abs(max(U)) < abs(min(V))) 
X a y’; 
else 
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y=X; 
end 
end 
count = 1; 
co = 1; 


while (gap(co) ~= 0) 
co=co+l1; 
count = count + 1; 


st = find(U == max(U)); 

if (length(st) > 1) % Randomly choose to break the tie 
i = ceil(rand(1)*length(st)); 
st = st(1); 

end 


x(st) = x(st) + 1; 
V=V+ S(st,:); 


nd = find(V == min(V)); 

if (length(nd) > 1) % Randomly choose to break the tie 
i= ceil(rand(1)*length(nd)); 
nd = nd(i); 

end 


y(nd) = y(nd) + 1; 


U=U + S¢,nd); 
JoTo To To Vo Vo To To To Vo Vo To Fo Fo Vo Vo 


xbarl = x(1:m)./sum(x(1:m)); 
ybar! = x(m+1:m+n)./sum(x(m+1:m-+n)); 
xbar2 = y(1:m)./sum(y(1:m)); 
ybar2 = y(m+1:m+n)./sum(y(m+1:m-+n)); 


lower(co) = max(min(xbar1'*A), min(xbar2*A)); 
upper(co) = min(max(A*ybar1), max(A*ybar2')); 
if (lower(co) < lower(co-1)) 

lower(co) = lower(co-1); 
end 
if (upper(co) > upper(co-1)) 

upper(co) = upper(co-1); 
end 


gap(co) = upper(co) - lower(co); 


if (mod(count,k) == 0) 
if (abs(max(U)) < abs(min(V))) 


V=-U; 
X=Y; 
else 
U=-V; 
y=x; 
end 
end 
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ah’ Marc trace EDR ; 
. Pe Coe eae 9 | 


if (gap(co) <= gapneed) 
break 
end 


% In order to save the memory space, 
% the following loop controls the result vectors (lower, upper and gap) 
% to have at most mx = 10,000 elements inside. 
if mod(count-1,mx) == 0 
r=r+1; 
disp([r, gap(co-1)}) %display the gap at every 10,000 iterations 
lower(1) = lower(co); 
lower = lower(1:mx); 
upper(1) = upper(co); 
upper = upper(1:mx); 
gap(1) = gap(co); 
gap = gap(1:mx), 
co= 1; 
end 
end 


compTime = cputime - t; 
treqflops = flops - f; 


if count > mx 
if co<5000 %number of iterations is in [r*10000 , r*10000+5000) 
lower = [lower((mx-5000+co+1):mx),lower(1:co)]; 
upper = [upper((mx-5000+co+1):mx),upper(1:co)]; 
gap = [gap((mx-5000+co+1):mx), gap(1:co)]; 
co = 5000; 
else Ynumber of iterations is in [r*10000+5000 , (r+1)*10000) 
lower = lower(1:co); 
upper = upper(1:co); 
gap = gap(1:co); 
end 
else %number of iterations < 10000 
lower = lower(1:co); 
upper = upper(i:co); 
gap = gap(1:co), 
end 


Upb = upper(end); 
Lwb = lower(end); 
gp = gap(end); 

P1 = xbar1; 

P2 = ybar2; 


return 
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APPENDIX F. MATLAB CODES FOR DMFP 


function [count, Upb,Lwb,gp,gap,compTime,reqflops,P1,P2] = DMFP(gapneed,A,w) 
% This function performs the new modified version of MFP algorithm in Gass, 

% Zafra's paper. 

r=0; 


k=1; 

{m,n] = size(A); 

alpha = 0.75; 

c = alpha*(max(max(A)) - min(min(A))); 

An=A+w; 

S = [zeros(m) An -c*ones(m, 1); 
-An’ zeros(n) c*ones(n,1); 
c*ones(1,m) -c*ones(1,n) 0]; 

[sm,sn] = size(S); 


U = zeros(sm,1); 
V = zeros(1,sn); 
x = zeros(sm, 1); 
y = zeros(1,sn); 


t = cputime; 
f = flops; 


%start the first iteration 

st = ceil(rand(1)*sm); % Row player randomly choose his row strategy 
x(st) = x(st) + 1; 

V=V+S(st,:); 


nd = find(V == min(V)); 

if (length(nd) > 1) % Randomly choose to break the tie 
1= ceil(rand(1)*length(nd)); 
nd = nd(i); 

end 


y(nd) = y(nd) + 1; 


U=U + S(:,nd); 
Jo To To To To To Vo To Mo To Yo Mo Yo Yo To Po 


xbarl = x(1:m)./sum(x(1:m)); % This part will cause zero division errors 
ybar! = x(m+1:m+n)./sum(x(m+1:m-+n)); % for some initial iterations, however MATLAB 
xbar2 = y(1:m)./sum(y(1:m)); % can continue to the end of mission. 


ybar2 = y(m+1:m+n)./sum(y(m+1:m-+n)); 


lower(1) = max(min(xbar1'*A), min(xbar2*A)); 
upper(1) = min(max(A*ybar1), max(A*ybar2’)); 
gap(1) = upper(1) - lower(1); 


if k == 
if (abs(max(U)) < abs(min(V))) 


X=y; 
else 
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U=-V'; 


y=X; 
end 
end 
count = 1; 
co=1: 


while (gap(co) ~= 0) 
co=co+l; 
count = count + 1; 


st = find(U == max(U)); 

if (length(st) > 1) % Randomly choose to break the tie 
i = ceil(rand(1)*length(st)); 
st = st(i); 

end 


x(st) = x(st) + 1; 
V=V+ S(st,:); 


nd = find(V == min(V)); 

if (length(nd) > 1) % Randomly choose to break the tie 
i = ceil(rand(1)*length(nd)); 
nd = nd(i); 

end 


y(nd) = y(nd) + 1; 


U=U + S¢,nd); 

Jo Vo To Fo To Jo To To Fo To Vo Yo To Mo Yo Yo 

xbar] = x(1:m)./sum(x(1:m)); % This part will cause zero division errors 

ybar! = x(m+1:m+n)./sum(x(m+1:m+n)); % for some initial iterations, however MATLAB 
xbar2 = y(1:m)./sum(y(1:m)); % can continue to the end of mission. 


ybar2 = y(m+1:m+n)./sum(y(m+1:m+n)); 
lower(co) = max(min(xbar1'*A), min(xbar2*A)); 
upper(co) = min(max(A*ybar1), max(A*ybar2’)); 
if (lower(co) < lower(co-1)) 

lower(co) = lower(co-1); 
end 
if (upper(co) > upper(co-1)) 

upper(co) = upper(co-1); 
end 
gap(co) = upper(co) - lower(co); 
if (mod(count,k) == 0) 

if (abs(max(U)) < abs(min(V))) 


V=-U'; 
x=y’; 
else 
U=-V;; 
y=X; 
end 
end 
if (gap(co) <= gapneed) 
break 
end 


% In order to save the memory space, 
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% the following loop controls the result vectors (lower, upper and gap) 
% to have at most mx = 10,000 elements inside. 
if mod(count-1,mx) == 0 
rer4+]; 
Z%oformat long e ; 
disp({r, gap(co-1)]) %display the gap at every 10,000 iterations 
lower(1) = lower(co); 
lower = lower(1:mx); 
upper(1) = upper(co); 
upper = upper(1:mx); 
gap(1) = gap(co); 
gap = gap(1:mx); 


co = 1; 
end 
% Dynamic scaling part 


JoTo To To To Yo To Vo To To Po To To To Po Lo Yo To Vo Vo Vo Vo 
if (mod(count,5000) == 0) % dynamically update w every 5000 iterations 
w = -0.5*(upper(co) + lower(co)); 
An=A+w; 
S = [zeros(m) An -c*ones(m, 1); 
-An’ zeros(n) c*ones(n,1); 
c*ones(1,m) -c*ones(1,n) 0]; 
end 
To To To To To Yo Vo To To To Vo To To To To Yo Mo To To Po Yo Mo Vo 
end 


compTime = cputime - t; 
reqflops = flops - f; 


if count > mx 
if co <5000 Y%number of iterations is in [r*10000 , r*10000+5000) 
lower = [lower((mx-5000+co+1):mx),lower(1:co)]; 
upper = [upper((mx-5000+co+1):mx),upper(1:co)]; 
gap = [gap((mx-5000+co+1):mx), gap(1:co)]; 
co = 5000; 
else %number of iterations 1s in [r*10000+5000 , (r+1)*10000) 
lower = lower(1:co); 
upper = upper(1:co); 
gap = gap(1:co); 
end 
else %number of iterations < 10000 
lower = lower(1:co); 
upper = upper(1:co); 
gap = gap(1:co); 
end 


Upb = upper(end); 
Lwb = lower(end); 
gp = gap(end); 


P1 = xbarl; 
P2 = ybar2; 
Return 
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